knitr::include_graphics("group.gif")Linear Mixed Model Analysis
1 INTRODUCTION
2 DATASET
dataset is taken from . the dataset is adopted from a research that objective is to study the relationship between time spent outdoors and mental wellbeing, across all of Scotland. The researcher contact all the Local Authority Areas (LAAs) and ask them to collect data for them, with participants completing the Warwick-Edinburgh Mental Wellbeing Scale (WEMWBS), a self-report measure of mental health and well-being, and being asked to estimate the average number of hours they spend outdoors each week.
2.1 Variables
the variables for the dataset are:
- local authority area (
laa): consist of 20 area - mental wellbeing score (
wellbeing) : Wellbeing score (Warwick Edinburgh Mental Wellbeing Scale). Range 15 - 75, with higher scores indicating better mental wellbeing - outdoor time (
outdoor_time) : Number of hours spent outdoors per week - density (
density) : Population density of local authority area (number of people per square km) - gender (
gender) : male or female - Participant Identifier (`ppt) : unique ID for each participants
individuals (ppt) are nested within district (laa). This means:
Level 1 is participants (ppt) Level 2 is district (laa)
3 Install and load packages
library(tidyverse)
library(ggplot2)
library(gtsummary)
library(readxl)
library(broom)
library(DT)
library(lme4)
library(kableExtra)
library(lmerTest)
library(lmtest)
library(broom.mixed)4 Read Data
data1 <- read_excel("D:/R Workspace/correlated numerical assignment/scotlandfinal.xlsx")
View(data1)5 Data wrangling
glimpse(data1)Rows: 132
Columns: 9
$ ppt <chr> "ID1", "ID10", "ID100", "ID101", "ID102", "ID103", "ID104…
$ name <chr> "Kendall Millar", "Dominik Hill", "Blair Harrison", "Zand…
$ laa <chr> "West Lothian", "Falkirk", "Falkirk", "Scottish Borders",…
$ outdoor_time <dbl> 20, 23, 29, 21, 10, 19, 13, 21, 16, 12, 22, 13, 26, 17, 2…
$ wellbeing <dbl> 37, 34, 39, 42, 37, 42, 38, 44, 47, 35, 42, 35, 29, 28, 2…
$ density <dbl> 421, 526, 526, 29, 19, 13, 31, 536, 480, 258, 44, 536, 42…
$ age <dbl> 74, 20, 77, 27, 68, 61, 73, 24, 29, 53, 82, 67, 30, 55, 5…
$ BMI <dbl> 22, 17, 28, 32, 26, 34, 19, 22, 36, 19, 31, 18, 19, 24, 2…
$ gender <chr> "Male", "Male", "Male", "Male", "Male", "Male", "Male", "…
summary(data1) ppt name laa outdoor_time
Length:132 Length:132 Length:132 Min. : 0.00
Class :character Class :character Class :character 1st Qu.:12.75
Mode :character Mode :character Mode :character Median :16.00
Mean :16.92
3rd Qu.:21.00
Max. :37.00
wellbeing density age BMI
Min. :16.00 Min. : 10.0 Min. :18.00 Min. :17.0
1st Qu.:34.00 1st Qu.: 19.0 1st Qu.:32.75 1st Qu.:21.0
Median :42.00 Median : 44.0 Median :55.50 Median :26.0
Mean :41.98 Mean : 477.5 Mean :54.01 Mean :26.3
3rd Qu.:48.00 3rd Qu.: 491.5 3rd Qu.:70.50 3rd Qu.:32.0
Max. :70.00 Max. :3581.0 Max. :85.00 Max. :36.0
gender
Length:132
Class :character
Mode :character
boxplot(data1$outdoor_time,
main = "Boxplot of Outdoor Time",
ylab = "Outdoor Time",
col = "lightblue",
border = "darkblue")data1 <- data1 %>%
mutate(gender = factor(gender, labels = c('male', 'female')))data1 <- data1 %>%
mutate(gender = factor(gender, labels = c('male', 'female')))data1<-data1 %>% mutate_if(is.character,~ as_factor(.))6 EDA
summarising the data
data1 %>%
select(laa, gender, outdoor_time, wellbeing, density) %>%
tbl_summary()| Characteristic | N = 1321 |
|---|---|
| laa | |
| West Lothian | 7 (5.3%) |
| Falkirk | 5 (3.8%) |
| Scottish Borders | 6 (4.5%) |
| Dumfries and Galloway | 6 (4.5%) |
| Argyll and Bute | 8 (6.1%) |
| Perth and Kinross | 7 (5.3%) |
| East Renfrewshire | 6 (4.5%) |
| Inverclyde | 5 (3.8%) |
| Midlothian | 6 (4.5%) |
| Moray | 7 (5.3%) |
| West Dunbartonshire | 6 (4.5%) |
| East Ayrshire | 7 (5.3%) |
| Shetland Islands | 7 (5.3%) |
| Stirling | 8 (6.1%) |
| City of Edinburgh | 8 (6.1%) |
| Orkney Islands | 6 (4.5%) |
| Na h-Eileanan Siar | 6 (4.5%) |
| Glasgow City | 8 (6.1%) |
| Highland | 7 (5.3%) |
| Angus | 6 (4.5%) |
| gender | |
| male | 19 (14%) |
| female | 113 (86%) |
| outdoor_time | 16 (13, 21) |
| wellbeing | 42 (34, 48) |
| density | 44 (19, 503) |
| 1 n (%); Median (Q1, Q3) | |
data1 %>%
ggplot(aes(x = outdoor_time, y = wellbeing)) +
geom_point() +
geom_smooth(method = lm)data1 %>%
ggplot(aes(x = outdoor_time, y = wellbeing,
col = gender, gender)) +
geom_point() +
geom_smooth(method = lm)7 Comparing groups with multilevel model
Start with null model (simplest model)
𝑠𝑐𝑜𝑟𝑒𝑖𝑗=𝛽0+𝑢0𝑗+𝑒𝑖𝑗
- 𝑠𝑐𝑜𝑟𝑒𝑖𝑗= the score attainment of participants 𝑖 in district𝑗
- 𝛽0 = the overall mean wellbeing score across district
- 𝑢0𝑗 = the effect of district𝑗on score. This is also level-2 residuals
- 𝑒𝑖𝑗 = individual-level residual. This is level-1 residuals
the null model : score
8 Single level analysis
using lm() as the outcome is assume as normally distributed
m.lm <- lm( wellbeing ~ 1, data = data1)
summary(m.lm)
Call:
lm(formula = wellbeing ~ 1, data = data1)
Residuals:
Min 1Q Median 3Q Max
-25.9848 -7.9848 0.0152 6.0152 28.0152
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 41.9848 0.9955 42.18 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.44 on 131 degrees of freedom
9 Multilevel analysis
We use the lme4 package, by starting with the constant only model (null model) with no explanatory variable or basically a random intercept with constant only model. the model is named as m0, the random effect is due to district.
m0 <-
lmer(wellbeing ~ 1 + (1 | laa),
data = data1, REML = FALSE)
summary(m0)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: wellbeing ~ 1 + (1 | laa)
Data: data1
AIC BIC logLik deviance df.resid
881.2 889.8 -437.6 875.2 129
Scaled residuals:
Min 1Q Median 3Q Max
-2.14222 -0.61901 0.02757 0.74655 1.95314
Random effects:
Groups Name Variance Std.Dev.
laa (Intercept) 99.92 9.996
Residual 27.24 5.219
Number of obs: 132, groups: laa, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 41.807 2.282 20.045 18.32 5.46e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
based on this model,the overall mean wellbeing score = 41.807. the mean for district j is estimated as 41.807 + 𝑈̂ 0𝑗 where 𝑈̂ 0𝑗 is the district residuals
the ICC is:
99.92 /(99.92 + 27.24)[1] 0.7857817
78.5% variance of the wellbeing score is attributed to the difference between district
or use tidy() for a proper table
tidy(m0) %>%
kbl() %>%
kable_styling()| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 41.807333 | 2.281684 | 18.32302 | 20.04486 | 0 |
| ran_pars | laa | sd__(Intercept) | 9.995826 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 5.218741 | NA | NA | NA | NA |
difference between multilevel and linear regression model
mlr <- lm(wellbeing ~ 1, data = data1)logLik(mlr) ; logLik(m0)'log Lik.' -508.4647 (df=2)
'log Lik.' -437.5814 (df=3)
the diiference is 70.88
anova(m0,mlr)the comparison is significant, thus the complex model of the the random intercept is the better model.
9.1 Random intercept model
9.1.1 adding the explanatory variable
We will model the effect of a individual-level variable outdoor time in the model
𝑠𝑐𝑜𝑟𝑒𝑖𝑗=𝛽0+𝛽1outdoor_time𝑖𝑗+𝑢0𝑗+𝑒𝑖𝑗
ri <- lmer(wellbeing ~ outdoor_time + (1 | laa),
data = data1,
REML = FALSE)
summary(ri)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: wellbeing ~ outdoor_time + (1 | laa)
Data: data1
AIC BIC logLik deviance df.resid
874.7 886.2 -433.4 866.7 128
Scaled residuals:
Min 1Q Median 3Q Max
-2.2347 -0.7223 0.1226 0.6432 1.8406
Random effects:
Groups Name Variance Std.Dev.
laa (Intercept) 100.71 10.036
Residual 25.24 5.024
Number of obs: 132, groups: laa, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 38.19177 2.59210 32.40431 14.734 6.44e-16 ***
outdoor_time 0.21340 0.07202 114.36148 2.963 0.00371 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
outdoor_tim -0.471
or
tidy(ri, conf.int = TRUE) %>%
kbl %>%
kable_styling()| effect | group | term | estimate | std.error | statistic | df | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 38.1917741 | 2.5920993 | 14.733916 | 32.40431 | 0.0000000 | 32.9144230 | 43.4691252 |
| fixed | NA | outdoor_time | 0.2133955 | 0.0720217 | 2.962931 | 114.36148 | 0.0037078 | 0.0707258 | 0.3560651 |
| ran_pars | laa | sd__(Intercept) | 10.0355904 | NA | NA | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 5.0235615 | NA | NA | NA | NA | NA | NA |
the equation for average fitted regression line (across district)
𝑠𝑐𝑜𝑟𝑒𝑖𝑗=38.1917741 + 0.2133955outdoor_time𝑖𝑗
the slope is fixed, while the intercept differs as it depends on the random effect of the district
9.1.2 Prediction
we can predict the wellbeing score based on the mixed model for each individual. the prediction is the average fitted regression plus the district intercept
pred_score <- fitted(ri)
head(pred_score, 10) 1 2 3 4 5 6 7 8
32.37244 31.74239 33.02276 40.05832 37.41243 43.90313 46.05396 44.41006
9 10
42.70001 33.07180
there will be 20 random effect because there were 20 district, and the random effect of each hospital as below
rand_ef <- ranef(ri)
head(rand_ef$laa, 20)to get the fitted values
ri_fitted <- augment(ri)ri_fitted %>%
slice(1:12)ri_fitted %>%
slice(42:52)confirmation with manual calculation
using the first observation:
- intercept: 38.1917741
- outdoor time : 20 minutes : 0.2133955(minutes)
- Community: West Lothian : -10.0872472
38.1917741-10.0872472 +(0.2133955*20)[1] 32.37244
the value from the table and manual calculation is similar which is 32.37244
9.1.3 Plot
ggplot(ri_fitted, aes(outdoor_time, .fitted, group = laa )) +
geom_point(alpha = 0.3) +
geom_line(alpha = 0.3) +
ylab('fitted wellbeing score') +
xlab('outdorr_time') +
ggtitle('The fitted value for random intercept model with outdoor time ') +
theme_bw()9.1.4 Variance
9.1.5 variance between district
in the constant only model the variance is 99.92 then the variance slightly inncrease after adding outdoor_time where model with outdoor_time as the explanatory variable now has the variance of 100.71 After accounting for coutdoor_time effects, the proportion of unexplained variance that is due to differences between district increases slightly to 100.71/(100.71+25.24)=80% .
100.71/(100.71+25.24)[1] 0.799603
9.1.6 within district variance
- constant only model variance is 27.24
- reduction of variance after adding outdoor_time
- model with outdoor_time as the explanatory variable 25.24
we can see that the addition of outdoor_time has increases the amount of variance at the district but not at the individual level. The between-district variance has increases from 99.92 to 100.71, and the within-school variance has reduced from 27.24 to 25.24.
10 Random slope model
now extend the random intercept model fitted before to allow both the intercept and the slope to vary randomly across district.
10.1 model
Wellbeing 𝑠𝑐𝑜𝑟𝑒𝑖𝑗=𝛽0+𝛽1outdoor_time𝑖𝑗+𝑢0𝑗+𝑢1𝑗outdoor_time𝑖𝑗+𝑒𝑖𝑗
rs <- lmer(wellbeing ~ outdoor_time + (1 + outdoor_time | laa),
data = data1, REML = FALSE)summary(rs)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: wellbeing ~ outdoor_time + (1 + outdoor_time | laa)
Data: data1
AIC BIC logLik deviance df.resid
866.6 883.9 -427.3 854.6 126
Scaled residuals:
Min 1Q Median 3Q Max
-2.26758 -0.59753 0.02406 0.62356 1.95808
Random effects:
Groups Name Variance Std.Dev. Corr
laa (Intercept) 44.2585 6.6527
outdoor_time 0.0891 0.2985 0.38
Residual 21.5661 4.6439
Number of obs: 132, groups: laa, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 38.27352 1.93169 19.04358 19.81 3.6e-14 ***
outdoor_time 0.21161 0.09662 13.78158 2.19 0.0462 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
outdoor_tim -0.231
rs <- lmer(wellbeing ~ outdoor_time + (1 + outdoor_time | laa), data = data1, control = lmerControl(optimizer = 'bobyqa'),
REML = FALSE)
summary(rs)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: wellbeing ~ outdoor_time + (1 + outdoor_time | laa)
Data: data1
Control: lmerControl(optimizer = "bobyqa")
AIC BIC logLik deviance df.resid
866.6 883.9 -427.3 854.6 126
Scaled residuals:
Min 1Q Median 3Q Max
-2.2676 -0.5977 0.0241 0.6236 1.9581
Random effects:
Groups Name Variance Std.Dev. Corr
laa (Intercept) 44.25665 6.6526
outdoor_time 0.08908 0.2985 0.38
Residual 21.56642 4.6440
Number of obs: 132, groups: laa, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 38.27359 1.93167 19.04249 19.81 3.6e-14 ***
outdoor_time 0.21161 0.09661 13.78262 2.19 0.0462 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
outdoor_tim -0.231
nicer output
tidy(rs) %>% kbl() %>%
kable_styling()| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 38.2735933 | 1.931666 | 19.813768 | 19.04249 | 0.0000000 |
| fixed | NA | outdoor_time | 0.2116052 | 0.096614 | 2.190214 | 13.78262 | 0.0462223 |
| ran_pars | laa | sd__(Intercept) | 6.6525669 | NA | NA | NA | NA |
| ran_pars | laa | cor__(Intercept).outdoor_time | 0.3800737 | NA | NA | NA | NA |
| ran_pars | laa | sd__outdoor_time | 0.2984660 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 4.6439655 | NA | NA | NA | NA |
10.2 The fitted values
rs_res <- augment(rs)
head(rs_res, 20)The outdoor effect for district𝑗is estimated as 0.2116051+𝑢̂ 𝑖𝑗, and the between-school variance in these slopes is estimated as 0.09.
That means for the average district we predict an increase of 0.2116051 points in the wellbeing score for every increase in 1 hour of outdoor time.
The intercept variance of 38.2735950 is interpreted as the between-district variance when outdoor_time=0
The intercept-slope correlation is estimated as 0.38 which means that district with a high intercept tend to have a steeper-than-average slope. This suggests that local authorities (laa) with higher baseline levels of wellbeing tend to show a stronger positive effect of outdoor_time on wellbeing. In other words, the more time people spend outdoors in districts that already have high average wellbeing, the greater the marginal gain in wellbeing from additional outdoor time.
10.3 Comparing models between random intercept and random slope
anova(ri, rs)There is very strong evidence that the outdoor time effect differs across district
10.4 Interpretation of random effects across district
The outdoor time effect for district 𝑗 is estimated as 0.21161+𝑈̂ 1𝑗, and the between-school variance in these slopes is estimated as 0.09.
For the average district we predict an increase of 0.21161 points in the wellbeing score for each hour of outdoor time. A 95% coverage interval for the school slopes is estimated as 0.21161±1.96√0.09 =−0.37639 to 0.79961.
Thus, assuming a normal distribution, we would expect the middle 95% of district to have a slope between −0.37639 and 0.79961.
ra.eff.rs <- ranef(rs, condVar = TRUE)
datatable(ra.eff.rs$laa)plot(ra.eff.rs)$laa
0 is equal to mean outdoor time of 16 hours per week
ra.eff.rs.sc <- ra.eff.rs$laa
names(ra.eff.rs.sc)[1] "(Intercept)" "outdoor_time"
ra.eff.rs.sc <- ra.eff.rs.sc %>%
rename(rs_slope = outdoor_time, rs_int = "(Intercept)")
ra.eff.rs.sc %>%
ggplot(aes( x = rs_int, y = rs_slope)) +
geom_point() +
geom_vline(xintercept = 0) +
geom_hline(yintercept = 0)This plot illustrates a positive correlation between intercepts and slopes:
districts with higher baseline wellbeing (positive random intercepts) tend to have stronger positive effects of outdoor_time on wellbeing (higher slopes).
Conversely, districts with lower-than-average baseline wellbeing (negative intercepts) tend to have weaker or negative relationships between outdoor_time and wellbeing.
#MODEL
expression(hat(wellbeingscore)[ij] ==
(38.2735950 + hat(u)[0 * j]) +
(0.21161 + hat(u)[1 * j]) * outdoor_time[ij])expression(hat(wellbeingscore)[ij] == (38.273595 + hat(u)[0 *
j]) + (0.21161 + hat(u)[1 * j]) * outdoor_time[ij])
plot(1, 1, type = "n", xlab = "", ylab = "", axes = FALSE)
text(1, 1,
expression(hat(wellbeingscore)[ij] ==
(38.2735950 + hat(u)[0 * j]) +
(0.21161 + hat(u)[1 * j]) * outdoor_time[ij]),
cex = 1.2)datatable(rs_res)manual calculation :
using the first observation:
- intercept: 38.2735950
- outdoor time : 20 minutes : 0.2116051(minutes)
- Community: West Lothian : -3.57859346508806
- outdoor time slope for west lothian :-0.321525576968577
38.2735950-3.57859346508806+((0.2116051-0.321525576968577)*20)[1] 32.49659
similar to the fitted value given in the table above: 32.49659260739616
#Adding a level 1 variable to the random slope model
rs_gen <- lmer(wellbeing ~ outdoor_time + gender + (1 + outdoor_time | laa),
data = data1, REML = FALSE,
lmerControl(optimizer = 'bobyqa'))
summary(rs_gen)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: wellbeing ~ outdoor_time + gender + (1 + outdoor_time | laa)
Data: data1
Control: lmerControl(optimizer = "bobyqa")
AIC BIC logLik deviance df.resid
867.7 887.9 -426.8 853.7 125
Scaled residuals:
Min 1Q Median 3Q Max
-2.21776 -0.61231 0.02179 0.64354 1.95652
Random effects:
Groups Name Variance Std.Dev. Corr
laa (Intercept) 43.74062 6.6137
outdoor_time 0.07373 0.2715 0.51
Residual 21.69867 4.6582
Number of obs: 132, groups: laa, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 37.17430 2.24156 32.86510 16.584 <2e-16 ***
outdoor_time 0.21329 0.09234 12.75828 2.310 0.0383 *
genderfemale 1.24578 1.26099 110.95959 0.988 0.3253
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) otdr_t
outdoor_tim -0.191
genderfemal -0.515 0.047
11 adding level2 explanatory variables
rs_den <- lmer(wellbeing ~ outdoor_time + density + (1 + outdoor_time | laa),
data = data1, REML = FALSE,
lmerControl(optimizer = 'bobyqa'))
summary(rs_den)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
method [lmerModLmerTest]
Formula: wellbeing ~ outdoor_time + density + (1 + outdoor_time | laa)
Data: data1
Control: lmerControl(optimizer = "bobyqa")
AIC BIC logLik deviance df.resid
867.7 887.9 -426.9 853.7 125
Scaled residuals:
Min 1Q Median 3Q Max
-2.26090 -0.56573 -0.00021 0.62588 2.03823
Random effects:
Groups Name Variance Std.Dev. Corr
laa (Intercept) 43.9730 6.6312
outdoor_time 0.0945 0.3074 0.24
Residual 21.4563 4.6321
Number of obs: 132, groups: laa, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 39.171950 2.157210 19.381769 18.159 1.25e-13 ***
outdoor_time 0.211476 0.098261 13.934109 2.152 0.0494 *
density -0.002048 0.002100 16.865581 -0.975 0.3433
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) otdr_t
outdoor_tim -0.283
density -0.446 0.028
anova(rs,rs_gen)anova(rs,rs_den)gender is not a significant covariates when included in the model.besides, model comparison also showed no significant difference between rs_den and rs_gen with the rs model.thus, the simpler model is chhosen. we remains with model with random slope for outdoor activity only
12 INTERACTION
as there is only one variable significant, no interaction checking needed
13 CHECKING FOR ASSUMPTION
Plot random effect
library(lattice)
randoms <- ranef(rs, condVar = TRUE)
dotplot(randoms)$laa
This plot above display the random effects estimates for different local authority areas “laa” in Scotland. The left side shows the variation in baseline levels of the outcome “wellbeing” score across different local authorities. Areas like Na h-Eileanan Siar and City of Edinburgh have higher-than-average baseline scores, while places such as Glasgow City and Falkirk fall below the overall mean, indicating lower baseline levels of the outcome.
The right panel, the “outdoor_time” illustrates the estimated random slopes for the effect of outdoor time across each local authority. Here, the estimates are clustered around zero, suggesting that the effect of outdoor time on the outcome is relatively uniform across regions or that the model has shrunk the estimates toward the overall fixed slope due to limited group-level variation. This indicates that while baseline outcomes vary substantially between areas, the influence of outdoor time on the outcome does not differ much by location.
plot(rs)the residual appear randomly scattered near the zero line .the assumption of homocedasticity was met
qqmath(rs)the points majority clustered on the line, so the residual is normally distributed. the assumption was met.
14 References
- Practical Linear Mixed Models. Kamarul Imran Musa. 25 March 20222.
- https://uoepsy.github.io/lmm/00_datasets.html